R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ──────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✔ ggplot2 3.3.5     ✔ purrr   0.3.4
✔ tibble  3.1.0     ✔ dplyr   1.0.9
✔ tidyr   1.2.0     ✔ stringr 1.4.0
✔ readr   2.1.2     ✔ forcats 0.5.1
package ‘tidyr’ was built under R version 4.0.5package ‘readr’ was built under R version 4.0.5── Conflicts ─────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ tidyr::expand() masks Matrix::expand()
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
✖ tidyr::pack()   masks Matrix::pack()
✖ tidyr::unpack() masks Matrix::unpack()

library(survival)
library(survminer)

Loading required package: ggpubr
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 methods overwritten by 'car':
  method                          from
  influence.merMod                lme4
  cooks.distance.influence.merMod lme4
  dfbeta.influence.merMod         lme4
  dfbetas.influence.merMod        lme4

library(cmprsk)

package ‘cmprsk’ was built under R version 4.0.5

library(tidyverse)
library(caret)

Loading required package: lattice

Attaching package: ‘caret’

The following object is masked from ‘package:survival’:

    cluster

The following object is masked from ‘package:purrr’:

    lift

library(survival)
library(survminer)
library(lubridate)

Attaching package: ‘lubridate’

The following objects are masked from ‘package:base’:

    date, intersect, setdiff, union

features = read_tsv("/Users/noa/Workspace/raxml_deep_learning_results/ready_raw_data/Pandit/ML/ML_edited_features.tsv")

New names:Rows: 825600 Columns: 129── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr   (7): msa_path, starting_tree_object, starting_tree_type, type, final_tree_topology, feature_msa_type, msa_name
dbl (121): ...1, Unnamed: 0, Unnamed: 0.1, index, Unnamed: 0.1.1, starting_tree_ind, spr_radius, spr_cutoff, final_ll, starti...
lgl   (1): starting_tree_bool
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

all_features<- read_tsv("/Users/noa/Workspace/raxml_deep_learning_results/ready_raw_data/Pandit/ML/all_features.tsv")

New names:Rows: 853120 Columns: 147── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr   (6): msa_path, starting_tree_object, starting_tree_type, type, final_tree_topology, feature_msa_type
dbl (141): ...1, Unnamed: 0, index, Unnamed: 0.1, starting_tree_ind, spr_radius, spr_cutoff, final_ll, starting_tree_ll, elap...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

all_features

all_features %>% distinct (feature_msa_pypythia_msa_difficulty, feature_mds_False_stress_20) %>% ggplot(aes(x=feature_msa_pypythia_msa_difficulty, y=feature_mds_False_stress_20))+geom_point()

all_features %>% distinct (feature_msa_pypythia_msa_difficulty, feature_mds_False_stress_3) %>% ggplot(aes(x=feature_msa_pypythia_msa_difficulty, y=feature_mds_False_stress_3))+geom_point()

Parsimony vs. Random

features %>% group_by(msa_path,starting_tree_ind,starting_tree_type) %>% summarise(is_global_max = mean(is_global_max)) %>% ggplot(aes(x= is_global_max, color =starting_tree_type))+geom_histogram()

`summarise()` has grouped output by 'msa_path', 'starting_tree_ind'. You can override using the `.groups` argument.

features %>% group_by(msa_path,starting_tree_ind,starting_tree_type) %>% summarise(mean_delta_ll = mean(delta_ll_from_overall_msa_best_topology)) %>% ggplot(aes(x= mean_delta_ll, fill =starting_tree_type))+geom_histogram()

`summarise()` has grouped output by 'msa_path', 'starting_tree_ind'. You can override using the `.groups` argument.

features %>% filter (starting_tree_type=='pars') %>% select (feature_mds_False_pca_0_3,feature_mds_False_pca_1_3,feature_mds_False_pca_2_3, delta_ll_from_overall_msa_best_topology, is_global_max) %>% ggplot(aes( x =feature_mds_False_pca_0_3,color=as.factor(is_global_max), ))+geom_density()

features %>% filter (starting_tree_type=='pars') %>% select (feature_mds_False_pca_0_3,feature_mds_False_pca_1_3,feature_mds_False_pca_2_3, delta_ll_from_overall_msa_best_topology, is_global_max) %>% ggplot(aes( x =feature_mds_False_pca_1_3,y = delta_ll_from_overall_msa_best_topology ))+geom_point()

NA

features %>% filter (starting_tree_type=='rand') %>% select (feature_mds_False_pca_0_3,feature_mds_False_pca_1_3,feature_mds_False_pca_2_3, delta_ll_from_overall_msa_best_topology, is_global_max) %>% ggplot(aes( x =feature_mds_False_pca_0_3,color=as.factor(is_global_max), ))+geom_density()

features %>% filter (starting_tree_type=='rand') %>% select (feature_mds_False_pca_0_3,feature_mds_False_pca_1_3,feature_mds_False_pca_2_3, delta_ll_from_overall_msa_best_topology, is_global_max) %>% ggplot(aes( x =feature_mds_False_pca_1_3,y = delta_ll_from_overall_msa_best_topology ))+geom_point()

features %>% filter (starting_tree_type=='rand') %>% select (feature_mds_False_pca_0_3,feature_mds_False_pca_1_3,feature_mds_False_pca_2_3, delta_ll_from_overall_msa_best_topology, is_global_max) %>% ggplot(aes( x =feature_mds_False_pca_1_3,y = delta_ll_from_overall_msa_best_topology ))+geom_point()

data = features %>% filter (starting_tree_type=='pars') %>% select (feature_mds_False_pca_0_3,feature_mds_False_pca_1_3,feature_mds_False_pca_2_3, delta_ll_from_overall_msa_best_topology, is_global_max,feature_msa_pypythia_msa_difficulty)
logistic_model <- glm(is_global_max ~ feature_msa_pypythia_msa_difficulty *feature_mds_False_pca_0_3*feature_mds_False_pca_1_3*feature_mds_False_pca_2_3, 
                      data = data, 
                      family = "binomial")
summary(logistic_model
        )

Call:
glm(formula = is_global_max ~ feature_msa_pypythia_msa_difficulty * 
    feature_mds_False_pca_0_3 * feature_mds_False_pca_1_3 * feature_mds_False_pca_2_3, 
    family = "binomial", data = data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9688  -0.9718  -0.4185   0.8960   2.6491  

Coefficients:
                                                                                                                  Estimate
(Intercept)                                                                                                        1.49377
feature_msa_pypythia_msa_difficulty                                                                               -7.48688
feature_mds_False_pca_0_3                                                                                         -0.42411
feature_mds_False_pca_1_3                                                                                          0.17615
feature_mds_False_pca_2_3                                                                                         -0.29570
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_0_3                                                      0.89976
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_1_3                                                     -0.57298
feature_mds_False_pca_0_3:feature_mds_False_pca_1_3                                                                0.11327
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_2_3                                                      0.60332
feature_mds_False_pca_0_3:feature_mds_False_pca_2_3                                                               -0.37910
feature_mds_False_pca_1_3:feature_mds_False_pca_2_3                                                               -2.83318
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_0_3:feature_mds_False_pca_1_3                           -0.50751
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_0_3:feature_mds_False_pca_2_3                            0.87724
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_1_3:feature_mds_False_pca_2_3                            7.16459
feature_mds_False_pca_0_3:feature_mds_False_pca_1_3:feature_mds_False_pca_2_3                                     -3.39414
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_0_3:feature_mds_False_pca_1_3:feature_mds_False_pca_2_3  8.69704
                                                                                                                  Std. Error
(Intercept)                                                                                                          0.04854
feature_msa_pypythia_msa_difficulty                                                                                  0.14110
feature_mds_False_pca_0_3                                                                                            0.05104
feature_mds_False_pca_1_3                                                                                            0.09061
feature_mds_False_pca_2_3                                                                                            0.08859
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_0_3                                                        0.14799
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_1_3                                                        0.24361
feature_mds_False_pca_0_3:feature_mds_False_pca_1_3                                                                  0.10303
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_2_3                                                        0.23684
feature_mds_False_pca_0_3:feature_mds_False_pca_2_3                                                                  0.10099
feature_mds_False_pca_1_3:feature_mds_False_pca_2_3                                                                  0.18815
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_0_3:feature_mds_False_pca_1_3                              0.27962
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_0_3:feature_mds_False_pca_2_3                              0.27288
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_1_3:feature_mds_False_pca_2_3                              0.47520
feature_mds_False_pca_0_3:feature_mds_False_pca_1_3:feature_mds_False_pca_2_3                                        0.22756
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_0_3:feature_mds_False_pca_1_3:feature_mds_False_pca_2_3    0.58465
                                                                                                                  z value
(Intercept)                                                                                                        30.774
feature_msa_pypythia_msa_difficulty                                                                               -53.063
feature_mds_False_pca_0_3                                                                                          -8.309
feature_mds_False_pca_1_3                                                                                           1.944
feature_mds_False_pca_2_3                                                                                          -3.338
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_0_3                                                       6.080
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_1_3                                                      -2.352
feature_mds_False_pca_0_3:feature_mds_False_pca_1_3                                                                 1.099
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_2_3                                                       2.547
feature_mds_False_pca_0_3:feature_mds_False_pca_2_3                                                                -3.754
feature_mds_False_pca_1_3:feature_mds_False_pca_2_3                                                               -15.058
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_0_3:feature_mds_False_pca_1_3                            -1.815
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_0_3:feature_mds_False_pca_2_3                             3.215
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_1_3:feature_mds_False_pca_2_3                            15.077
feature_mds_False_pca_0_3:feature_mds_False_pca_1_3:feature_mds_False_pca_2_3                                     -14.915
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_0_3:feature_mds_False_pca_1_3:feature_mds_False_pca_2_3  14.876
                                                                                                                  Pr(>|z|)    
(Intercept)                                                                                                        < 2e-16 ***
feature_msa_pypythia_msa_difficulty                                                                                < 2e-16 ***
feature_mds_False_pca_0_3                                                                                          < 2e-16 ***
feature_mds_False_pca_1_3                                                                                         0.051905 .  
feature_mds_False_pca_2_3                                                                                         0.000844 ***
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_0_3                                                      1.2e-09 ***
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_1_3                                                     0.018671 *  
feature_mds_False_pca_0_3:feature_mds_False_pca_1_3                                                               0.271620    
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_2_3                                                     0.010852 *  
feature_mds_False_pca_0_3:feature_mds_False_pca_2_3                                                               0.000174 ***
feature_mds_False_pca_1_3:feature_mds_False_pca_2_3                                                                < 2e-16 ***
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_0_3:feature_mds_False_pca_1_3                           0.069528 .  
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_0_3:feature_mds_False_pca_2_3                           0.001306 ** 
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_1_3:feature_mds_False_pca_2_3                            < 2e-16 ***
feature_mds_False_pca_0_3:feature_mds_False_pca_1_3:feature_mds_False_pca_2_3                                      < 2e-16 ***
feature_msa_pypythia_msa_difficulty:feature_mds_False_pca_0_3:feature_mds_False_pca_1_3:feature_mds_False_pca_2_3  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 536277  on 412799  degrees of freedom
Residual deviance: 430414  on 412784  degrees of freedom
AIC: 430446

Number of Fisher Scoring iterations: 4

features %>% distinct (msa_path, feature_msa_pypythia_msa_difficulty)

all_features %>% distinct(msa_path, feature_msa_pypythia_msa_difficulty)

features %>% distinct (msa_path, feature_msa_pypythia_msa_difficulty)

all_features %>% group_by(msa_path,starting_tree_ind,starting_tree_type,feature_mds_False_pca_1_3,feature_mds_False_pca_2_3) %>%
filter(msa_path=='/groups/pupko/noaeker/data/New_MSAs/Pandit_msas/PF08033.fasta') 

NA

features %>% distinct (msa_path, feature_msa_pypythia_msa_difficulty)

features %>% group_by(msa_path,starting_tree_ind,starting_tree_type,feature_mds_True_pca_0_3,feature_mds_True_pca_1_3,feature_mds_True_pca_2_3) %>% summarise(is_global_max = max(delta_ll_from_overall_msa_best_topology<=0.1)) %>% filter (starting_tree_type=='rand') %>% filter(msa_path=='/groups/pupko/noaeker/data/New_MSAs/Pandit_msas/PF08021.fasta')  %>% ggplot(aes(x=feature_mds_True_pca_1_3, y = feature_mds_True_pca_2_3, color =as.factor(is_global_max), shape = starting_tree_type))+geom_point()

`summarise()` has grouped output by 'msa_path', 'starting_tree_ind', 'starting_tree_type', 'feature_mds_True_pca_0_3', 'feature_mds_True_pca_1_3'. You can override using the `.groups` argument.

features  %>% mutate(is_global_max = max(delta_ll_from_overall_msa_best_topology<=0.1)) %>% ggplot(aes(x=running_time, y = feature_mds_True_pca_2_3, color =as.factor(is_global_max), shape = starting_tree_type))+geom_point()

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.